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\  . 

Dynamic  stall  effects  are  analyzed  in  this  investigation  for  cases  of 
Ci(4->-an  inertially  static  airfoil  in  a  flow  field  rotating  at  constant  rate 
(gust  response),  and  (2)  an  airfoil  pitching  at  constant  rate  in  a  steady 
flow  field.  The  method  used  is  a  boundary  layer  solution  of  the  momentum- 
integral  equation  by  a  modified  von  Karman-Pohlhausen  technique. 

Work  accomplished  by  Docken  in  1982  using  this  method  to  match 
Kramer's  experimental  results  for  gust  response  is  reviewed,  corrected, 
and  continued.  The  validity  of  the  closure  equation  and  the  assumptions 
key  to  its  derivation  are  examined,  concluding  that  the  closure  equation 
is  justified.  A  better  match  of  Kramer's  airfoil  sections  results  in 
dynamic  stall  predictions  very  close  to  experimental  data.  The  effect 
of  varying  airfoil  thickness  and  camber  is  investigated. 

By  consideration  of  the  non-Newtonian  motion  of  the  boundary  layer 
on  the  surface  of  a  pitching  airfoil,  the  momentum-integral  method  is 
extended  to  the  second  case.  Using  the  Moore-Rott-Sears  model  for  flow 
separation  criteria,  analytical  results  were  computed  and  compared  with 
experimental  data.  Reduction  in  adversity  of  the  pressure  gradient 
accounts  for  only  a  fraction  of  the  total  dynamic  effect,  and  it  is 
proposed  that  mass  introduction  into  the  boundary  layer  from  the  free- 
3 cream  may  be  a  strongly  contributing  factor.  This  phenomena  is 
demonstrated  to  have  a  large  effect,  and  an  argument  is  presented  for 
the  proper  amount  of  mass  introduction. 


INVESTIGATION  OF  EFFECTS  CONTRIBUTING  TO  DYNAMIC  STALL 


USING  A  MOMENTUM- INTEGRAL  METHOD 

I.  Introduction 

Discussion 

Dynamic  stall  is  a  phenomenon  that  has  received  considerable  attention 
over  the  past  two  decades.  It  is  understandably  of  intense  interest  to  the 
helicopter  and  compressor  industries,  where  much  of  the  research  has  been 
conducted,  and  more  recently  has  received  the  attention  of  the  aerospace 
industry  in  applications  to  winged  space  shuttle  vehicles  (Ref  23:2). 

The  majority  of  the  research  has  been  experimental,  and  most  of  the 

analytical  work  has  relied  on  some  sort  of  experimental  data  as  input. 

Further,  because  of  the  importance  of  practical  application  (i.e.,  helicopter 
blades),  nearly  all  of  the  work  in  recent  years  has  dealt  with  oscillatory 
motion  of  the  airfoil.  Consequently,  there  is  still  much  to  be  learned 
and  understood  about  the  basic  physical  flow  phenomena  which  contribute  to 
dynamic  stall. 

The  concept  is  not  difficult.  When  an  airfoil  angle  of  attack  is 
increased  rapidly  through  a  range  that  includes  the  static  stall  angle, 
the  maximum  lift  and  drag  are  significantly  increased,  stall  is  delayed  to 

a  higher  angle  of  attack,  and  the  corresponding  loss  of  lift  is  often  much 

more  severe  in  nature.  The  problem  then  is  one  of  complexity.  Contributing 
factors  are  both  viscous  and  inviscid.  The  phenomenon  is  strongly  dependent 


on  the  rate  of  pitch,  free-stream  velocity,  and  frequency  and  amplitude 
(for  an  oscillating  airfoil),  and  shows  varying  degrees  of  dependence  on 
airfoil  shape,  3D  flow  effects,  Reynolds  number,  and  Mach  number 
(Ref  1:304). 

Problem  Statement 

It  is  far  beyond  the  scope  of  this  study  to  attempt  to  analytically 
predict  dynamic  stall  effects  on  an  airfoil.  Rather,  the  purpose  of  this 
thesis  is  to  attain  a  better  physical  understanding  of  some  of  the  u.tri- 
buting  causes  of  dynamic  stall.  As  such,  restriction  is  made  to  t 
analysis  of  boundary  layer  effects  for  an  airfoil  undergoing  a  cor  rate 

of  change  in  angle  of  attack.  It  would  seem  essential  that  the  phyai.cs  of 
this  relatively  simple  problem  should  be  well  understood  and  documented 
prior  to  an  analytical  study  of  more  complex  oscillatory  motions,  but  such 
has  not  been  the  case. 

The  method  of  analysis  will  be  to  integrate  a  solution  of  the  unsteady 
(time-dependent)  boundary  layer  equations,  using  a  modified  von  Karman-Pohl- 
hausen  momentum-integral  technique,  until  stall  is  indicated,  and  determine 
the  maximum  angle  of  attack  at  this  point.  Stall  of  the  airfoil  is  here 
defined  as  the  condition  of  flow  separation  at  quarter-chord  of  the  airfoil, 
with  flow  separation  being  initially  defined  as  the  condition  of  zero  shear 
at  the  surface  as  indicated  by  the  velocity  profile  shape  parameter  (this 
definition  of  separation  is  later  modified  using  the  MRS  conditions  for  the 
case  of  a  pitching  wing).  To  use  this  technique,  the  classical  boundary 
layer  assumptions  must  be  applied. 

Two  distinct  cases  will  be  analyzed.  The  first  is  that  of  an  inertially 


static  airfoil  in  a  rotating  flow  field,  and  this  will  heretofore  be 
referred  to  as  the  "gust-response"  problem.  The  second  is  that  of  a 
pitching  airfoil  in  a  steady  flow  field,  and  will  be  referred  to  as  the 


"pitching  airfoil"  problem.  These  may  appear  at  first  to  be  the  same 
problem,  differing  only  in  the  reference  frame  of  the  observer.  Experi¬ 
mental  results,  however,  indicate  greatly  differing  effects,  as  will  be 
shown.  The  analysis  of  this  report  will  also  point  out  the  differences 
between  these  two  cases. 

Background 

Max  Kramer  reported  the  results  of  an  experimental  study  of  the  gust- 
response  problem  in  1932.  By  fixing  a  wing  section  between  a  set  of 
moveable  guide  vanes  and  a  diffuser,  he  was  able  to  take  lift  measurements 
while  constantly  varying  the  angle  (with  respect  to  chordline)  of  the 
free-stream  velocity.  Based  on  his  results,  Kramer  established  a  linear 
relationship  between  the  increase  of  maximum  lift  coefficient  (ovei  the 
steady  flow  case)  and  the  non-dimensional  pitch  rate  (Ref  2). 

An  analytic  study  of  the  gust-response  problem  conducted  by  Richard 
Docken  in  1982  successfully  predicted  a  relationship  in  close  agreement 
with  Kramer  (Ref  3).  He  developed  a  modified  von  Karman-Pohlhausen  integral 
method  for  unsteady  flow  and  applied  it  to  the  boundary  layer  of  a  represen¬ 
tative  Joukowski  airfoil.  The  outstanding  feature  of  this  method  is  the 
development  of  a  closure  equation  to  handle  the  unknown  unsteady  term. 

Two  experimental  studies  of  the  pitching  airfoil  problem  are  of 
particular  interest.  In  1978  Arthur  Deekens  and  William  Kuebler  conducted 
smoke  tunnel  tests  of  an  airfoil  section  pitching  at  various  constant  rates 


and  m  various  free-stream  steady  velocities.  Motion  pictures  were  taken 
of  the  flow,  allowing  observation  of  flow  separation  at  quarter-chord  versus 
instantaneous  angle  of  attack.  By  comparison  with  static  conditions,  they 


were  able  to  determine  a  direct  correlation  between  the  delay  in  stall 


Deekens /Kuebler 


ts  of  Experimental  Investigations  of  Dynamic  Stall 
Kramer,  Deekens/Kuebler ,  and  Daley 


II .  Gust  Response  Problem 

Review  of  Docken's  Work 

It  is  notable  that  in  Docken's  approach  to  the  problem  of  gust  response 
each  step  of  this  theory  development  and  development  of  the  method  of 
integration  is  simply  an  unsteady  modification  of  the  corresponding  steady 
flow  model  (Ref  6:17,20).  This  idea  is  reinforced  by  the  fact  that,  when 
his  program  is  applied  to  a  case  of  zero  pitch  rate,  time  derivatives 
become  zero,  and  all  equations  and  the  resulting  solution  reduce  to  the 
steady-state  case. 

The  theory  development  leading  to  the  unsteady  boundary  layer 
momentum-integral  equation  is  reviewed  in  Appendix  A.  In  short,  the 
principles  of  continuity  and  momentum  are  applied  to  an  incremental  control 
volume  in  the  boundary  layer.  The  unsteady  Euler's  equation  is  then 
employed  to  substitute  for  the  pressure  gradient  term.  By  careful  manipu¬ 
lation  of  terms,  the  momentum-integral  equation  for  unsteady  flow  results: 
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Docken's  modification  to  the  von  Karman-Pohlhausen  integral  method  to 


handle  unsteady  flow  is  reviewed  in  Appendix  B.  As  in  the  well  known 
steady  flow  model  (Ref  7:158-160),  the  displacement  thickness  and  momentum 
thickness  are  substituted  into  the  momentum-integral  equation.  The  velocity 


profile  is  defined  as  a  polynomial  expression  of  v/6,  and,  by  applying 
known  boundary  conditions,  the  velocity  profile  can  be  related  to  a  single 
shape  parameter: 
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Again  as  in  the  steady-flow  model,  additional  parameters  are  defined 
and  substituted  into  the  governing  equation,  and  the  "working"  equations 
for  stepwise  integration  become 
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The  single  exception  to  this  parallelism  of  the  unsteady  flow  derivation 
to  the  steady  flow  derivation  arises  from  the  introduction  of  the  displace¬ 
ment  and  momentum  thicknesses  into  the  momentum-integral  equation.  A 
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transient  term  of  the  form  — —  developed,  requiring  a  closure  equation. 

O  t 

The  necessary  equation  was  derived  as 


(Ref  3) 
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The  limits  of  applicability  of  the  assumptions  associated  with  this 
derivation  were  not  investigated  in  depth,  however,  and  both  the  derivation 
and  the  underlying  assumptions  will  be  examined  in  Section  III  of  this 
report . 

The  method  used  to  step-wise  numerically  integrate  Eqs  (3)  and  (4) 
along  the  surface  of  a  Joukowski  airfoil  is  detailed  in  Appendix  C. 
Coordinates  on  the  airfoil  and  velocity  at  the  outer  edge  of  the  boundary 
layer  are  obtained  by  a  Joukowski  transformation  from  flow  about  a  circular 
cylinder,  the  unsteadiness  being  approximated  as  pseudo-steady  Joukowski 
solutions  at  the  changing  angles  of  attack.  Necessary  derivatives  at  each 
step  are  computed  using  standard  difference  methods,  and  the  shape  parameter 
can  be  examined  at  any  given  point  on  the  airfoil  to  determine  a  point  of 
boundary  layer  separation  (Ref  6:23-25). 

Docken  defined  aerodynamic  stall  as  the  condition  of  flow  separation 
at  the  quarter-chord  of  the  airfoil.  Flow  separation  was  determined  to 
occur  when  shear  at  the  airfoil  surface  equals  zero,  as  indicated  by  the 
shape  parameter  (A  »  -12,  or  K  *  -.1567).  Using  this  criteria,  the  above 
method  of  solution  was  applied  to  a  J015  airfoil.  By  assuming  the 
C]~vs.-a  slope  remains  linear  and  equal  to  2ir/radian,  and  applying  three- 
dimensional  flow  corrections,  Docken 's  analytical  results  were  expressed  as 


AC. 


.343 


Ca 


comparing  remarkably  well  with  Kramer's  experimental  results: 
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The  significance  of  this  result  is  not  only  that  an  approximate  integral 
method  had  been  modified  to  give  such  a  degree  of  accuracy,  but  that,  like 
the  experimental  results,  the  relationship  between  gust  pitch  rate  and 
dynamic  stall  angle  of  attack  is  linear  (Ref  6:27).  Also,  and  just  as 
important,  the  cause  of  the  increased  angle  of  attack  can  be  attributed  to 
a  distinct  physical  effect  —  a  reduction  in  the  adversity  of  the  pressure 
gradient  as  determined  from  Euler's  equation. 

Modification  and  Results 

Upon  close  examination  of  the  Docken  program,  an  error  was  discovered 
in  the  Joukowski  transformation.  The  error  was,  in  fact,  causing  the  method 
of  solution  to  be  applied  to  an  airfoil  with  thickness  less  than  the  desired 
15%.  This  error  was  corrected  and  some  other  minor  modifications  were  made 
to  the  original  work. 

Upon  running  the  solution  again  for  a  J015  airfoil,  the  data  remained 
linear,  but  with  a  different  linear  coefficient.  The  corrected  data  is 
displayed  graphically  in  Figure  2. 

Converting  this  data  in  order  to  compare  with  Kramer's  results,  the 


results  of  Figure  2  may  be  expressed  as 
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In  radians  then,  and  correcting  the  form  of  the  non-dimensional  pitch 


rate, 
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Assuming,  by  classical  airfoil  theory,  a  2ir/radian  curve  slope  (Ref  8:125), 
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Since  Kramer's  experiment  dealt  with  a  wing  with  aspect  ratio  approxi¬ 
mately  equal  to  five,  a  three-dimensional  flow  correction  is  required.  By 
wing-section  theory,  the  lift  curve  slope  per  degree  of  the  corresponding 
wing  may  be  obtained  by 


a 


a 


(Ref  9:11) 
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III.  Closure  Equation 


Review  of  Derivation 

Prior  to  Docken's  work,  the  von  Karman-Pohlhausen  technique  had  not 
been  used  to  solve  an  unsteady  flow  problem,  because  the  unsteady  terms 
in  the  momentum-integral  equation  required  an  additional  unknown  closure 
equation.  The  derivation  of  the  required  equation  is  simple  and  concise, 
and  the  assumptions  appear  to  be  well-founded,  at  least  superficially 
(Ref  6:14-15).  In  order  to  investigate  the  validity  of  these  assumptions, 
and  hence  the  equation,  and  to  then  further  examine  the  impact  on  the 
problem  solution,  it  is  necessary  first  to  fully  understand  the  derivation 

To  find  the  equation  of  closure,  it  is  first  observed  that,  for  any 
location  and  time  in  the  flow,  the  displacement  thickness  is  related  to 
the  boundary  layer  thickness  as 
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Also,  for  laminar  flows,  the  boundary  layer  thickness  is  related  to 
the  velocity  of  the  potential  flow  as 
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Then,  taking  the  time  derivatives  of  Eqs  (6)  and  (7), 


3$,  a*  3C 

3t  C1  3t  F  6  (8) 


and 


35 

3t 


1/2  au 

It 


+ 


(=^) 
U  x 
e 


1/2 


(9) 


Substituting  Eqs  (7)  and  (9)  into  Eq  (8)  gives 


36 

3? 


1  »  - 


C  C 

-LI  (JL) 

2U  x' 

e  e 


1/2  3U  3C0 

— S.  +  c  _ -  (  v  ) 

3t  1  3t  U  x' 
e 


1/2 


3C1  V  1/2 

*  c2  sr  <uV 

e 


(10) 


Now,  the  restrictive  assumption  is  made  that  the  time  dependence 
of  both  and  is  small  compared  to  that  of  U  ,  and  Eq  (10)  becomes 
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Applying  Eqs  (6)  and  (7)  to  Eq  (11)  then  yields  the  required  closure 


equation: 
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Investigation  of  Assumptions 

3C1  3C2 

The  key  assumption  then  is  that  and  are  negligibly  small, 
and  the  question  arises  as  to  what  magnitudes  are  acceptable  to  the 
solution.  From  Eqs  (10)  and  (11),  it  can  be  seen  that  the  closure 
equation  is  valid  if 


C  C 

-LI  (LL) 

2U  U  x ' 
e  e 


1/2  3U 


3t 


_e  a 


1/2  3U 


SC, 


.-1.  2  (.y.  )  '  ~  _le  +  c  III  (_v  \ 

2U  U  x  3t  1  3t  U  x 

e  e  e 


1/2 


SC. 


+  r  _ -  (  —  ) 

Z  St  x 
e 


1/2 
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Then,  to  compare  the  relative  magnitudes  of  the  time  derivatives, 


Eq  (13)  requires  that 


From  this  it  can  be  seen  that,  while  this  relationship  is  certainly 


satisfied  if  the  original  assumption  is  true,  it  is  also  satisfied  if  the 
3C2  3C1  . 

sum,  ,  is  comparatively  small.  This  could  certainly  be 

the  case  if  the  time  variations  of  and  are  of  opposite  sign. 

Either  of  these  assumptions  could  be  evaluated  by  the  degree  to  which 
Eq  (14)  is  satisfied.  The  velocity  and  acceleration  terms  are  available  at 
any  given  point  on  an  airfoil  by  the  Pohlhausen  method  previously  discussed; 
however,  a  scheme  must  be  developed  to  determine  C^,  C^,  and  the  corresponding 
time  derivatives . 

Consider  first  the  term  C^.  By  definition,  as  in  Eq  (6),  is  the 
ratio  of  displacement  thickness  to  boundary  layer  thickness.  This  can  be 
written  in  terms  of  the  Pohlhausen  velocity  profile  shape  parameter  as 


(15) 


Then  the  time  derivative  of  ,  also  in  terms  of  the  shape  parameter, 


is 


1  3A 

120  3t 


(16) 


When  an  airfoil  is  at  an  angle  of  attack  such  that  separation  is 
indicated  at  the  quarter-chord  (i.e.,  the  case  under  consideration),  the 


shape  parameter  varies  from  +7.052  at  the  stagnation  point  to  -12.0  at 

quarter-chord  (Ref  7:211);  hence,  for  a  static  airfoil  and  —  for  a 

dx  Dt 

pitching  airfoil  are  negative  quantities.  The  fact  that  a  pitching  airfoil 

demonstrates  delayed  separation  would  infer  that  the  associated  ^  term  is 

<3t 

causing  A  to  be  less  negative,  thus  is  a  positive  term.  If  this  is 

3C 

indeed  the  case,  then  7^—  is  a  negative  term,  by  Eq  (16). 

Since  the  shape  parameter  is  available  in  the  solution  at  each  point 

on  the  airfoil,  C  is  easily  determined.  The  term  ,  and  consequently 
9C1  . 

—  ,  can  be  calculated  in  the  following  manner.  The  total  change  in  the 
O  t 

shape  parameter  can  be  calculated  between  two  given  points  on  the  airfoil 

by  both  the  steady  solution  (time  independent)  and  the  pitching  solution 

3A 

(time  dependent).  For  the  steady  case,  -^  =  0 .  Defining  the  total 
derivative  as 


3A  3s  3A 

3t  +  3t  3s 


let  (— )  be  approximated  by  a  steady,  or  time-independent,  difference  of 

A  between  two  given  points  in  close  proximity,  and  (— )  be  approximated 

C  2 

by  the  time-dependent  difference  between  the  same  two  points.  Then 


, pA '  .DA-.  /  3A  ,  3s  3A\  .  3s  3A\  _  3A 

''Dt  ^  "  ut '  x  V3t  3t  3s;  ”  3t  3s  "  3t 
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Consider  now  the  term  C^.  If  the  velocity  of  the  potential  flow  over 
a  wedge  is  proportional  to  a  power  of  the  length  coordinate  as 


U(x)  =  kxm 


(18) 


where  m  is  a  specific  function  of  the  included  wedge  angle,  then  the 
solutions  of  Falkner  and  Skan  show  that  fcr  any  given  wedge  there  is  a 
constant  such  that 


5 


s  (uV 

e 


1/2 


(Ref  7:164-166) 


(19) 


The  constant  is  only  a  function  of  the  geometry;  i.e.,  wedge  angle. 


Comparing  Eqs  (7)  and  (19),  it  can  be  seen  that  if  such  self-similar 
solutions  can  be  applied  to  an  airfoil  approximated  by  a  series  of  panels 
(i.e.,  the  wing  is  made  up  of  a  series  of  short  linear  ramps),  each  making  a 


specific  "wedge  angle"  with  the  free-stream,  then  C  could  be  determined. 


3C_ 


Further,  since  C  is  only  a  function  of  geometry,  tt ~  on  a  pitching 

L  at 

airfoil  could  also  be  determined  through  the  change  in  slope  of  the  panel 
relative  to  the  free-stream  over  an  increment  of  time. 


Accuracy  of  Closure  Equation 

Evaluation  of  the  closure  equation  assumptions  was  accomplished  through 


8 


the  determination  of  C  ,  C  ,  tt—  ,  and  - -  on  a  typical  pitching  airfoil 

1  2  o  t  o  t 

as  it  approached  separation  at  the  quarter-chord.  Using  the  unsteady 

3U 

Pohlhausen  method,  U^,  ,  and  A  are  readily  available.  An  independent 

set  of  Falkner-Skan  solutions  was  computed,  and  a  curve  fit  of  the  results 
made,  such  that  the  constant  could  be  determined  for  any  included 
wedge  angle.  This  and  the  techniques  described  in  the  previous  paragraphs 
were  then  incorporated  into  the  computer  solution,  such  that  at  each  point 
on  the  airfoil  a  comparison  could  be  made  of  all  of  the  above  terms. 
Further,  each  of  the  following  relationships  was  computed,  and  a  comparison 
made : 


C  C  3U 
12  e 


(20a) 


C  C  3U  3C„ 

1  2  e  .  p  2  ,  p 

~  5T  C1JT*C! 

e 


(20b) 


From  Eq  (13)  it  can  be  seen  that  the  degree  to  which  the  computed 
value  of  (20a)  approaches  that  of  (20b)  is  representative  of  the  accuracy 
of  the  closure  equation. 

Investigation  was  made  of  eight  separate  cases,  each  varying  in 
free-stream  velocity  and  pitch  rate,  such  that  non-dimensional  pitch  rates 
varied  from  .009  to  .035.  The  following  observations  were  made: 


3C  3U 

1  9  6 

a.  t —  and  t —  are  indeed  smaller  than  — ,  but  not 

'  .  .  9C1C  9S 

negligibly  so.  The  magnitudes  of  t — =-  and  r —  increase 


3U 

with  increasing  chord  position,  while  ■  is  initially 

increasing,  then  decreasing.  The  ratio  of 
3U  SC,  3C 

rr- -  :  t— : —  :  —  varied  from  104:102:1  to  30:15:1,  the 

3t  3t  3t 

latter  (worst)  case  occurring  at  the  quarter-chord 
(separation  point)  for  the  situation  of  slow  free-stream 
velocity  and  slow  pitch  rate. 


30  .  3C2 

b.  —  was  negative  in  value,  while  C  ,  C  ,  and  ttt-  were 

8v..1_  K 


3t 


positive.  The  products  and  C2  —  then  had  a 

canceling  effect  on  one  another. 


& 


c.  The  comparison  of  terms  (20a)  and  (20b)  is  shown  in 
Figure  3  and  was  nearly  identical  for  all  test  cases. 
This  would  indicate  an  accuracy  of  the  closure  equation 
as  reflected  in  Figure  4. 


It  may  be  noted  here  that  the  validity  of  the  method  used  to  determine 

9S 

^2  3n<*  "3t~  *-8  1uest:'-on>  since  it  is  not  known  if  the  velocity  on  the 
airfoil  taken  from  Joukowski  transformation  agrees  with  Eq  (18).  To 
investigate  this,  a  program  was  written  to  simultaneously  and  independently 
compute  two  velocities  on  a  representative  pitched  airfoil.  One  velocity 
was  that  obtained  from  the  Joukowski  transformation,  while  the  other  began 
with  an  identical  value  at  the  leading  edge  of  the  airfoil  and  progressed 
over  a  series  of  incremental  panels  (representative  of  wedges  of  identical 
slope  to  the  airfoil  as  discussed  earlier).  This  was  done  using  Eq  (18), 
where  x  is  the  additive  distance  traversed  over  the  panels  from  the 


20 


Chord 


stagnation  point  and  m  is  determined  from  the  slope  of  the  panel  relative 
to  the  free-stream.  The  relationship  is 


m 


(Ref  7:155) 


where  Yir  is  the  included  angle  of  the  wedge.  The  constant  k  can  be  determined 

for  each  panel  at  the  initial  juncture  by  letting  U  at  the  end  of  one  panel 

e 

equal  at  the  beginning  of  the  next.  Thus,  for  a  juncture  of  panels  i 
and  i  +  1  at  position  x 


The  results  of  this  study  are  shown  in  Figure  5.  The  "panel"  velocity 
varied  in  error  from  6%  (negative)  near  the  leading  edge  to  15%  (positive) 
at  the  quarter-chord.  The  velocities  were  nearly  identical  in  the 
vicinity  of  .05  chord. 

Consider  the  hypothetical  effect  of  forcing  the  panel-derived  velocities 
to  agree  with  the  Joukowski-derived  velocities.  This  would  require  "wedge" 
panels  with  an  increased  included  angle  near  the  leading  edge,  and  a 
decreased  included  angle  aft  of  about  .05  chord.  The  effect  would  be 


3C 

to  reduce  the  value  of  C„  and  r —  forward  of  .05  chord,  and  increase  the 

2  ot 


values  aft  of  .05  chord. 


Joukowski  transform 


Figure  5.  Comparison  of  Potential  Velocities  on  the  Surface  of  an  Airfoil 


Now  consider  the  effect  on  the  results  of  Figure  3.  It  can  be  shown 

3C2  . 

that  r —  is  affected  considerably  more  than  C  by  any  change  in  the 

at  2 

included  angle.  Then  the  end  effect  would  reduce  the  relative  error  of 
the  closure  equation  everywhere  on  the  airfoil.  Note  also  that  Figure  3 
indicates  that  the  closure  equation,  like  the  above  velocities,  has  near 
zero  error  in  the  vicinity  of  .05  chord. 

In  order  to  determine  the  effect  that  error  in  the  closure  equation 
has  on  the  prediction  of  dynamic  change  in  stall  angle  of  attack,  errors 
of  50%  (very  conservative,  as  shown  in  Figure  3)  and  25%  (more  realistic  but 
still  conservative  as  argued  above)  were  introduced  into  Eq  (12).  The 
unsteady  Pohlhausen  equations  were  then  derived  in  the  same  manner  used 
by  Docken.  The  results  are  that  a  50%  error  in  the  closure  equation 
results  in  a  13.54%  under-prediction  of  the  dynamic  effect  on  stall  angle 
of  attack,  and  25%  error  in  the  closure  equation  results  in  under-predicting 
this  effect  by  7.04%.  The  error  presented  then  in  the  closure  equation 
has  kept  the  final  results  on  the  conservative  side  and  shows  that  the 
closure  equation  is  a  good  one,  but  for  slightly  different  reasons  than 
stated  by  Docken. 


Further  Investigation  of  Gust  Response 


Better  Match  of  Kramer  s  Airfoils 


In  a  closer  examination  of  Kramer  s  experimental  work,  it  may  be  noted 
that  data  was  taken  from  two  considerably  different  airfoil  profiles  —  a 
Gottingen  459  (symmetric)  airfoil  and  a  Gottingen  398  (cambered)  airfoil. 
The  resulting  dynamic  effects  on  maximum  coefficient  of  lift  were  almost 
identical,  and  Kramer  states  that  the  profile  of  the  wing  apparently  has 
little  influence  (Ref  2). 

The  nature  of  the  unsteady  solution  method  developed  by  Docken  readily 
lends  itself  to  an  investigation  of  this  hypothesis.  To  this  end,  a  simple 
computer  program  was  written  to  determine  what  specific  input  data  were 
necessary  to  reproduce  desired  variations  in  thickness  and  camber  to  a 
Joukowski  airfoil.  It  was  found  that  a  12.64%  thick  symmetrical  airfoil 
very  closely  matched  the  Gottingen  459  airfoil  in  shape,  and  a  Joukowski 
airfoil  with  13.72%  thickness  and  3.17%  camber  likewise  closely  matched 
the  Gottingen  398  airfoil.  A  comparison  of  airfoil  shapes  is  shown  in 
Figure  6  and  7. 

After  a  minor  modification  to  the  Docken  program  to  allow  for  camber 
variation,  computer  solutions  were  obtained  for  maximum  angle  of  attack 
of  each  of  the  above  airfoils.  Nine  sets  of  conditions  were  run  for  each 
airfoil,  and  the  results  are  shown  in  Figure  8. 

The  results,  when  converted  as  in  Section  II  to  compare  with  Kramer's 


data,  are  as  follows: 


For  the  Gottingen  459  (symmetric)  airfoil: 


es  Used 


Gottingen  459 


mge  in  Separation  Angle  of  Attack 
Response  for  Gottingen  Airfoils 


1.283  —  corrected  for  R.  =  5 


For  the  Gottingen  398  (cambered)  airfoil: 

f-337  r 


.241  —  corrected  for  R  =  5 


For  comparison,  recall  that  Kramer's  experimental  results  for  each 


of  the  above  airfoil  profiles  were 


And  for  Docken's  J015  airfoil: 

f-3"  r 


.220—  corrected  for  R 


Note  that  with  the  better  Joukowski  shape  approximation  of  the 
Gottingen  airfoils,  results  were  considerably  closer  to  Kramer's  solution. 
It  is  meaningless  to  discuss  relative  error  between  any  of  the  computer 
results  and  Kramer's  experimental  results  because  of  the  assumptions 
inherent  to  the  momentum-integral  method  of  solution  and  approximations 
necessary  to  compare  results;  however,  the  proximity  of  the  above  solutions 


is  certainly  remarkable.  Note  also  that  the  profile  of  the  airfoil  did, 
in  fact,  affect  the  results. 

Thickness  and  Camber  Effects 

Kramer  discusses  measurements  only  on  wings  with  profiles  of  the 
two  airfoils  examined  here,  and  apparently  bases  his  assumption  that  profile 
has  little  effect  only  on  those  results.  From  Figure  8  it  can  be  seen  that 
these  particular  airfoils  do  not  differ  much  in  change  of  stall  angle  ^f 
attack,  even  though  the  airfoil  shapes  are  considerably  different.  The 
fact  that  there  does  exist  a  difference  based  on  shape  alone,  however, 
suggests  further  investigation  is  in  order. 

Two  studies  were  conducted  regarding  airfoil  shape.  First,  using  a 
symmetrical  Joukowski  airfoil,  the  thickness  was  varied,  with  a  linear 
coefficient  of  increase  of  stall  angle  of  attack  per  non-dimensional  pitch 
rate  being  calculated  for  maximum  thicknesses  of  10%,  12.68%,  15%,  20%, 
and  25%.  Secondly,  a  15%  thick  Joukowski  airfoil  was  varied  in  camber, 
and  the  same  coefficient  calculated  for  maximum  cambers  of  0%,  1%,  2.5%, 

5%,  and  7.5%.  The  results  are  displayed  graphically  in  Figures  9  and  10, 
respectively. 

The  dynamic  stall  effect  was  decreased  as  airfoil  thickness  increased 
and  as  camber  increased.  While  camber  has  a  relatively  slight  effect, 
the  thickness  effect  was  dramatic.  It  would  appear  then  that  airfoil 
profile  is  indeed  an  important  factor  in  dynamic  stall  effect.  It  should 
be  noted  that,  while  Kramer's  airfoils  appear  very  different,  they  differ 
by  only  about  1%  in  maximum  thickness. 


Effect  of  Airfoil  Camber  on  Gust  Response  Solution 


V.  Pitching  Airfoil  Problem 


Non-Inertial  Control  Volume  Analysis 


As  discussed  previously,  the  only  readily  apparent  difference  between 
the  gust-response  problem  and  the  pitching  airfoil  problem  is  the  frame 
of  reference.  It  was  anticipated  that,  due  to  the  similarity  of  these 
problems,  it  may  be  possible  to  investigate  the  dynamic  stall  effects  on 
a  pitching  airfoil  by  paralleling  the  approach  taken  in  the  gust-response 
solution.  The  procedure  then  is  to  derive  a  momentum- integral  equation 
for  the  boundary  layer  by  developing  the  continuity  and  momentum  equations 
for  a  fluid  element  in  the  boundary  layer. 

The  continuity  equation  presents  no  problems.  In  exactly  the  same 
analysis  as  for  the  gust-response  solution,  the  principle  of  continuity 
is  applied  to  an  incremental  control  volume,  with  the  resulting  equation 
(see  Appendix  A): 


=  '  h  /  Dudydx 


In  order  to  use  the  momentum  equation  on  this  control  volume ,  however 
certain  modifications  must  be  made.  The  momentum  equation,  derived  from 
Newton' 8  Second  Law,  is  usually  applied  to  a  control  volume  that  is  either 
fixed,  or  translating  at  constant  velocity,  with  respect  to  an  inertial 
reference  frame.  The  surface  of  the  pitching  airfoil,  a.id  thus  the 
attached  control  volume  in  the  boundary  layer,  is  neither.  By  an  analysis 


of  particle  dynamics,  certain  "hypothetical  body  forces"  can  be  determined 
which  account  for  the  non-Newtonian  motion  of  the  control  volume,  and  with 
the  addition  of  these  forces  the  control  volume  analysis  may  proceed  as 
if  it  were  fixed  in  inertial  space  (Ref  10:109-114).  The  derivation  and 
examination  of  the  hypothetical  body  forces  is  included  in  Appendix  D. 
These  body  forces  in  the  x-direction  are: 


dB  *  -  dm  |  +  82v  -  3^*1 

where : 

*  Ra2  cos  <t>  , 

B2  =  2a  , 


and  R  is  the  length  from  the  point  of  rotation  of  the  airfoil  to  a  reference 
point  in  the  control  volume,  and  4>  is  the  angle  between  line  R  and  the 
local  surface  tangent  of  the  airfoil. 

Now  referring  to  Figure  11,  the  sum  of  the  forces  in  the  x-direction 


can  be  equated  to  the  net  rate  of  transport  of  momentum  out  of  the  control 
volume  plus  the  flux  of  momentum  in  the  control  volume,  such  that: 


Figure  11.  Momentum  Conservation  in  a  Non-Inertial 
Boundary  Layer  Control  Volume 
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(22) 


Note  that  the  classical  boundary  layer  assumption  is  made  that  a 
pressure  gradient  exists  only  in  the  x-direction.  If  Lq  (21)  is  now 
substituted  into  (22),  and  dividing  by  pdx,  the  result  is 


o 


3(uu) 

3x 


dy  -  U 

e 


(23) 


Consider  the  hypothetical  body  force  terms  in  Eq  (23).  B^dx  is  very 
small  in  comparison  to  B^;  i.e., 


a2  dx  «  R  a2  cos  $ 
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It  can  be  shown  that  R  cos  <j>  is  a  large  term  except  very  near  the 
leading  edge.  On  the  airfoil  under  consideration,  R  cos  <J>  >  1  for  chord 
positions  between  .006  and  .250.  Since  dx  «  1 ,  it  may  be  assumed  that 


B^dx) 


dy 


It  is  demonstrated  in  Appendix  D  that  a  potential  flow  velocity  close 
.  .  * 

to  and  tangent  to  the  airfoil  surface,  and  of  magnitude  U  as  observed 
in  the  inertial  frame  of  reference,  will  appear  to  an  observer  fixed  in 
the  control  volume  as 


i  ^  •  A  •  A 

U  ■  (U  -  R  a  sin  <J>)  X  -  (r  a  cos  <J>)  y 


(24) 


Note  that  there  is  a  constant  y-component  of  velocity  toward  the  surface 
of  the  airfoil;  but  if  the  surface  is  impermeable,  then  this  component  is 
necessarily  zero  at  y  *  0.  Since  this  difference  must  be  compensated  for 
in  the  boundary  layer,  it  is  assumed  that  this  y-component  of  velocity 
is  a  linear  function  for  0  ^  y  <  S  .  Then 


R  ct  cos  $ 


V 


-  R  a  cos  <j>  , 


y  >  <5 


Now, 
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h  h 

f  dy  ~  82  f  vdy 

J  ft 
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R  a  cos  cj) 


h 

J  ydy  +  J  (-  R  a  cos  (J>)dy 


Carrying  out  the  integrations ,  this  reduces  to 


(h  -  6)  R  a2 


cos  4> 


The  total  effect  of  the  hypothetical  body  forces  on  a  control  volume 
where  the  upper  limit  is  the  edge  of  the  boundary  layer  (h  =  6)  is  zero. 
Although  at  this  point  h  is  still  arbitrary,  the  boundary  conditions  will 
eventually  be  applied  at  y  *  6,  therefore  these  terms  may  be  dropped  from 


i 


Eq  (23).  In  the  interest  of  argument,  the  effect  of  these  terms 
not  disappear  if  v  does  not  vary  linearly  in  the  boundary  layer. 


does 

Cases 
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of  parabolic  and  elliptic  variations  were  investigated,  with  the  resulting 
(small)  terms  being  carried  completely  through  the  Pohlhausen  derivation 
and  then  applied  to  an  airfoil.  The  effects  on  the  final  solution  were 
negligible  in  both  cases. 

Eq  (23)  thus  becomes: 


3(uu) 


d?  -  ue  /  I; dy  *  h  /  udy  * 


-i  r"fi  dy-^ 

P  J  dx  p 


At  this  point  the  unsteady  EuleT's  equation  is  used  to  substitute 
for  the  unknown  pressure  gradient.  Letting  an  asterisk  denote  velocity 
as  viewed  from  the  inertial  reference  frame, 


From  Eq  (24),  an  observer  in  the  control  volume  will  see  a  magnitude 


U  -  U  -Rasin<j> 


Therefore,  in  terns  of  the  velocity  seen  by  the  control  volume, 
Eq  (26)  becomes 
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Substituting  equation  (27)  for  the  pressure  gradient  in  Eq  (25) 


yields 
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With  the  exception  of  the  (added)  second  term  on  the  rij,ht-hand  side 


of  this  equation,  it  is  exactly  the  same  as  that  derived  for  the  gust- 
response  problem,  and  these  terms  may  be  manipulated  as  shown  in 
Appendix  A.  The  resulting  modified  momentum-integral  equation  is 
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Modification  of  von  Karman-Pohlhausen  Method 

Again,  as  with  the  gust-response  problem  in  Appendix  B,  this  can  be 
written  in  terms  of  the  displacement  and  momentum  thicknesses.  If  the 
integral  of  the  first  term  is  evaluated  with  the  outer  edge  of  the  boundary 
layer  as  the  upper  limit,  the  result  is 
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Multiplying  this  equation  by  — —  gives 


(R  a  sm  $)  —  - — 

U  3x 
e 


+ 


U  —  — • 

e  V  dx 


3U 

<5  <5  2  . 

3U 

6,  6  2  36, 

e 

12  1 

e 

1  2  1 

3x 

<57  v  u 

2  e 

3t 

<$2  v  3t 

T  6 
w 


2 


UU 

e 


Applying  the  closure  equation  (12),  and  rearranging  terms,  this  may  be 
written  as 
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From  this  point,  the  derivation  of  the  working  equations  for  the 
unsteady  Pohlhausen  method  closely  resembles  that  for  the  gust-response 
problem  as  outlined  in  Appendix  B;  therefore,  only  the  differences  will 
be  highlighted  here. 

Due  to  the  change  in  appearance  of  the  unsteady  pressure  term,  the 


corresponding  boundary  condition  becomes 
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The  resulting  working  equations  for  stepwise  integration  are: 
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Results  for  Pitching  Motion 


The  existing  program  for  the  gust-response  solution  was  modified 
such  that  the  non-Newtonian  motion  of  the  airfoil  could  be  accounted  for. 
Eqs  Ola)  and  (31b)  were  implemented,  which  necessitated  computation  of 
the  new  terms  f^(K)  and  R  a  sin  <J> .  In  addition,  recall  that  the  velocity 
U  and  its  corresponding  derivatives  must  now  be  those  that  are  observed 


in  a  frame  of  reference  fixed  to  the  airfoil  surface.  A  subroutine  to 
the  initial  program  was  thus  written  to  compute  R  a  sin  4>  (geometry 
dependent)  and  to  correct  the  Joukowski  transformation-derived  velocity 
for  motion  of  the  airfoil. 

Upon  computing  the  change  in  maximum  angle  of  attack  for  various 
combinations  of  free-stream  velocity  and  pitch  rate  and  using  the 
separation  condition  of  zero  shear  stress  at  the  wall,  it  was  found 
that  the  results  remained  linear  with  respect  to  non-dimensional  pitch 
rate.  The  results  are  shown  in  Figure  13,  being  the  set  of  data  points 
labeled  "pitching  motion".  Note  that  the  effect  of  a  pitching  airfoil 
is  approximately  three  times  greater  than  that  of  an  inertially-f ixed 
airfoil  in  a  pitching  airflow  (gust  response).  The  cause  of  this  increase 
in  angle  of  attack  of  separation  appears  to  be  due  to  the  effectively 
reduced  velocity  of  the  free-stream  by  varying  amounts,  depending  on  the 
location  on  the  airfoil.  This  velocity  modification  both  delays  the 
onset  of  and  relaxes  the  severity  of  the  adverse  pressure  gradient. 

Moore-Rott-Sears  Model  and  Results 

In  considering  the  pitching  airfoil  as  a  "moving  wall"  type  of  problem, 
it  is  also  necessary  to  reevaluate  the  criteria  for  determining  the  point 
of  separation.  Up  to  this  point,  separation  has  been  assumed  to  occur 
when  shear  at  the  wall  vanishes,  which  for  a  static  wall  is  the  limit 
between  forward  and  reverse  flows  in  the  boundary  layer.  For  a  wall 
moving  in  the  direction  of  flow,  this  is  not  true  (i.e.,  the  flow  may 
obtain  a  velocity  less  than  the  velocity  of  the  wall  but  greater  than  zero). 
The  MRS  (Moore-Rott-Sears)  model  provides  criteria  of  separation  based 
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on  the  velocity  profile  which  will  account  for  motion  of  the  wall 
(Ref  11:123-126).  It  defines  separation  as  the  condition  where  both 
velocity  and  shear  become  zero  in  a  singular  fashion  at  some  point  in 
the  boundary  layer,  as  seen  by  an  observer  in  the  frame  of  reference  of 
the  wall.  This  condition  is  illustrated  in  Figure  12. 

A  separate  program  was  written  to  provide  a  table  of  values  relating 
a  separation  A  (or  K)  for  any  given  value  of  Then,  by  computing 

^wall^e  3t  <luarter“c*lor<*  anc*  comparing  with  this  table,  the  dynamic-stall 
condition  could  be  determined  by  iteration  until  the  final  A  matches  the 
indicated  tabular  A.  It  should  be  noted  that  this  is  only  valid  for 
A  >  -17.76,  since  the  shape  parameter  K  reaches  a  mathematical  minimum 
at  that  point  (Ref  12:295-296). 

Application  of  this  new  stall  criteria  had  a  considerable  effect 
on  the  solution.  The  results  are  displayed  in  Figure  13,  labeled 
"MRS  model".  Note  that  the  dynamic  effect  on  angle  of  attack  is  more  than 
twice  that  using  the  old  criteria,  and  the  total  dynamic  effect  is  about 
seven  times  that  of  the  gust-response  solution.  It  is  also  of  interest 
that  the  change  in  angle  of  attack  is  now  no  longer  linear  with  non-dimen¬ 
sional  pitch  rate,  but  is  curved  similar  to  the  data  in  Figure  1;  however, 
the  extent  of  the  dynamic  effect  falls  well  short  of  the  experimental  data. 

Mass  Introduction  and  Results 

The  combined  dynamic  effects  thus  far  discussed  predict  an  angle  of 
attack  change  considerably  less  than  experimental  results  indicate,  as  shown 
in  Figure  13.  This  implies  that  other  contributing  factors  are  involved, 
and  it  was  suggested  that  mass  introduction  into  the  boundary  layer  may 
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undary  Layer  Velocity  Profiles  At  and  Near  Separation 


be  a  contributing  cause.  Note  that  Eq  (24)  indicates  that  a  potential 
flow  velocity  at  the  edge  of  the  boundary  layer  does  in  fact  have  a 
component  toward  the  surface  when  viewed  from  the  control  volume.  While 
this  normal  component  is  everywhere  small,  a  review  of  boundary  layer 
solutions  involving  "suction"  and  "blowing"  indicates  that  very  small  mass 
flow  rates  can  have  dramatic  effects  on  boundary  layer  separation 
(Ref  7:380-399). 

To  implement  this  idea  into  the  method  of  solution,  Eq  (21)  was 
amended  to  account  for  an  additional  increment  of  mass,  C(R  a  cos  cj> )  dx, 
into  the  top  of  the  control  volume  (where  C  is  some  as  yet  undetermined 
constant) : 


m 

top 


pudydx  - 


C(R  a  cos  4>)  dx 


(32) 


The  substitution  of  this  into  the  momentum  equation,  with  all  other 


terms  carried  forward  as  previously  described,  results  in  an  amended 


momentum-integral  equation: 
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Casting  Eq  (33)  into  the  Pohlhausen  form  as  before  presented  little 
problem,  and  all  of  the  shape  parameters  and  related  functions  remained 
the  same.  The  additional  term  altered  the  working  equation  (31a)  only. 
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Two  problems  now  become  apparent.  First,  some  value  must  be  assigned 
to  the  constant  C.  Secondly,  the  appearance  of  in  the  new  term  is 
undesirable,  and  renders  Eq  (34)  unusable  in  its  present  form. 

Consider  the  fact  that,  by  injecting  mass  into  the  top  of  the  boundary 
layer,  a  stagnation  flow  problem  is  being  imposed  upon  the  boundary  layer. 
The  impinging  velocity  of  a  two-dimensional  potential  flow  near  the 
stagnation  point  is  given  as 


■—  =  -  ay  (Ref  7:96) 


Applied  to  the  problem  at  hand,  this  would  mean  that  at  the  outer  edge 


of  the  boundary  layer  the  velocity  desired  is: 

V  3  -  a<S  (R  a  cos  <J>) 

e 

The  constant  a  is  dependent  upon  the  dimensions  of  the  problem.  In 
Appendix  E  it  is  argued  that,  for  the  Joukowski  airfoil  problem,  a  reasonable 
guess  for  the  value  of  a  is  four.  Using  this, 


C  =  46 


and  Eq  (34)  becomes 
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Using  Eqs  (31b)  and  (36)  as  the  working  equations,  solutions  were 
again  computed  for  various  free-stream  velocities  and  pitch  rates.  The 
results  are  shown  in  Figure  13,  labeled  "Mass  Introduction  (46)".  The 
injection  of  this  small  amount  of  mass  into  the  boundary  layer  has  more 
than  doubled  the  dynamic  increase  in  stall  angle  of  attack.  This  is 
attributed  to  the  fact  that  the  injected  mass  is  accelerating,  or  imparting 
additional  energy  to,  the  fluid  in  the  boundary  layer,  hence  retarding 
separation. 

The  computer  program  used  to  calculate  these  solutions  is  included 
for  reference  in  Appendix  F. 

The  justification  for  the  use  of  a  =  4  in  equation  (35)  that  is 
outlined  in  Appendix  E  is  somewhat  arbitrary,  and  certainly  arguments 
could  be  presented  for  other  values  of  this  constant.  The  possibility 
was  considered  that  careful  selection  of  this  value  might  provide  data 
that  would  match  experimental  data.  In  this  interest,  various  values 
were  tried.  Figure  13  shows  the  results  for  a  =  20,  labeled  "Mass 
Introduction  (206)".  It  is  apparent  that  the  experimental  results  cannot 
be  matched  in  this  manner.  It  is  interesting  and  possibly  significant  to 
note  that,  for  the  selected  value  a  =  4,  the  slope  of  the  plotted  data 
becomes  closely  tangent  to  the  slope  of  the  experimental  data. 

This  implies  that  still  other  factors,  in  addition  to  mass  introduc¬ 
tion,  might  be  contributing  to  the  dynamic  stall  effect,  although  it  is 
clearly  implied  that  mass  introduction  is  an  important  contributor. 


Effects  Contributing  to  Change  in  Separation  Angle  of  Attack 
Due  to  the  Pitching  Motion  of  an  Airfoil 


VI.  Conclusions  and  Recommendations 


Conclusions 

For  the  case  of  an  airfoil  encountering  a  gust  of  constant  o  ,  the 
increase  in  dynamic  stall  angle  of  attack  can  be  attributed  almost 
exclusively  to  a  reduction  in  the  adversity  of  the  pressure  gradient  on 
the  upper  surface  of  the  airfoil. 

The  use  of 
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is  acceptable  for  use  as  an  equation  of  closure  for  the  unsteady  von 
Karman-Pohlhausen  integral  method.  A  conservative  estimate  of  the  error 
inflicted  by  its  use  causes  an  under-prediction  of  increase  in  dynamic 
stall  angle  of  attack  by  less  than  10%,  and  the  analysis  seems  to  indicate 
that  the  approximation  is  really  much  better  than  this  estimate. 

The  increase  in  dynamic  stall  angle  of  attack  for  gust  response  is 
inversely  proportional  to,  and  a  strong  function  of,  the  maximum  thickness 
of  the  airfoil.  It  is  inversely  proportional  to,  and  a  weak  function  of, 
the  maximum  camber  of  the  airfoil. 

Contributing  factors  to  the  increase  in  dynamic  stall  angle  of  attack 
for  an  airfoil  pitching  at  constant  a  in  a  steady  flow  are: 


a.  A  reduction  in  the  adversity  of  the  pressure  gradient 
on  the  upper  surface  of  the  airfoil. 


b.  An  effectively  reduced  free-stream  velocity  due  to 
the  downstream  motion  of  the  surface  of  the  airfoil, 
causing  further  relaxation  of  the  severity  of  the 
adverse  pressure  gradient. 

c.  A  delay  in  boundary  layer  separation  due  to  the  fact 
that  the  motion  of  the  airfoil  surface  allows  the 
velocity  gradient  to  build  past  the  condition  of  zero 
wall  shear  to  the  condition  of  zero  velocity  and  shear 
at  some  point  in  the  boundary  layer. 

d.  An  acceleration  of  the  fluid  in  the  boundary  layer  due 
to  mass  introduction  into  the  top  of  the  boundary  layer. 

The  factors  above  may  not  be  all  of  the  contributing  effects,  indicat¬ 
ing  that  other  phenomena  should  be  considered  to  account  for  the  total 
dynamic  effect.  However,  it  is  clear  that  the  mass  introduction  is  an 
important  contributor  to  the  effect. 

Recommendations 

In  considering  the  idea  of  mass  introduction  into  the  boundary  layer 
of  a  pitching  airfoil,  the  selection  of  46  as  the  appropriate  constant 
was  a  rough  estimate  at  best,  and  intended  only  to  show  trend  on  the 
solution.  An  in-depth  study  of  the  concept  of  mass  introduction  is 
appropriate . 

In  this  report,  the  effect  of  the  shedding  of  starting  vortices  due 
to  the  change  in  the  airflow  pattern  about  the  airfoil  has  been  assumed  to 
have  negligible  contribution  to  the  solution.  Work  currently  in  progress 
by  K.  Tupper  in  this  area  would  indicate  that  this  effect  should  be 
considered  (Ref  14).  A  study  combining  the  work  of  this  report  with 
that  of  Tupper  seems  a  good  step  toward  identifying  other  possible  dynamic 


stall  contributors. 
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Appendix  A:  Derivation  of  the  Unsteady  Momentum-Integral  Equation 


The  momentum-integral  equation  for  unsteady  flow  is  derived  by 
integrating  the  equation  of  motion  over  the  thickness  of  the  boundary 
layer.  The  equation  of  motion  in  turn  is  derived  by  applying  the 
principles  of  continuity  and  conservation  of  momentum  in  the  direction 
of  flow  to  an  incremental  control  volume  in  the  boundary  layer. 

The  principle  of  continuity  (mass  conservation)  states  that  the 
net  rate  of  mass  flow  across  the  control  surface  must  equal  the  time 
rate  of  mass  reduction  within  the  control  volume.  Referring  to  Figure  14, 
it  is  shown  that 

'  h  3  h  -  h 

ratop  +[o/  pudy  +  3^  /  pudydxr  J  eudy 

~  h 

■  '  37  /  pd? 

o 

For  incompressible  flow,  i.e.  p  =  constant,  the  last  term  in  this 
expression  is  zero.  Therefore 


3  r 

-  J  pudydx 


o 


(Al) 
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.ass 


The  principle  of  momentum  conservation  (in  the  x-direction)  requires 


that  the  sum  of  the  forces  acting  on  the  fluid  in  the  control  volume  is 
equal  to  the  net  rate  of  transport  of  momentum  across  the  control  surface 
plus  the  time  flux  of  momentum  within  the  control  volume.  Referring  to 
Figure  15,  it  is  shown  that 
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Simplifying  and  rearranging  terms, 
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Substituting  Eq  (Al)  for  the  unknown  term  in  (A2), 


Note  that,  since  the  integrals  are  taken  in  the  y-direction,  the 
partial  derivatives  with  respect  to  dimension  x  (and  also  time)  can  be 
arbitrarily  brought  inside  the  integrals.  In  addition,  by  boundary  layer 
assumptions,  the  pressure  gradient  in  the  x-direction  is  assumed  to  be 
the  total  (spatial)  pressure  gradient.  Employing  these  concepts,  and 
dividing  the  entire  equation  by  pdx,  Eq  (A3)  becomes 


/  b  (uu)dy  '  °« 


The  Euler's  equation  of  motion  for  unsteady  flow  can  be  shown  to  be 


1  dP 


U 
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(Ref  27) 


This  may  be  applied  to  Eq  (A4)  in  substitution  for  the  pressure 
gradient.  After  rearranging  terms,  Eq  (A4)  now  becomes 
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Consider  the  time  variant  terms  of  Eq  (A5): 
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Further, 


since  U  is 
e 


independent  of  y, 
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3t 
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Now  consider  the  spatially  variant  terms  of  Eq  (A5): 


(A5) 


( A6) 
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Substitution  of  Eqs  (A6)  and  (A7)  into  Eq  (A5)  yields  the  momentum- 


integral  equation  for  unsteady  flow: 
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Appendix  B:  Development  of  the  von  Karman-Pohlhausen 
Method  for  Unsteady  Flow 


By  comparing  the  unsteady  momentum-integral  equation  (A8)  with  the 
definitions  of  the  boundary  layer  displacement  and  momentum  thickness  in 
their  integral  form  (Ref  7:140-141),  it  is  evident  that,  for  arbitrary  h, 
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Expanding  the  first  and  third  terms  and  rearranging  yields 
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At  this  point  an  equation  of  closure,  Eq  (5),  is  used  to  substitute 
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for  the  difficult  term  7^—  .  This  closure  equation  and  its  validity  is 
discussed  at  length  in  Section  III.  The  result  is 
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Finally,  by  adding  and  subtracting  the  quantity 
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the  momentum-integral  equation  becomes 
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The  velocity  profile  of  the  boundary  layer  must  satisfy  the  following 


boundary  conditions: 


at  y  •  0, 

u  =  0 

(B2a 

a2,. 

3U  3u 

a  u 

V  x-r 
3y 

=  -  (U  r-2-  +  r^) 

e  dx  dt 

(B2b 

at  y  *  6,  u  =  U 
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(B2c 


(B2d 


(B2e 


Since  the  existence  of  five  boundary  conditions  will  allow  for  the 
solution  of  five  free  constants,  the  velocity  profile  may  be  expressed  as 
a  fourth-degree  polynomial,  where  H 

-  a  +  Bn  +  cn2  +  Dn3  +  En4  (B3) 

e 

Applying  (B2a)  yields 

A  -  0  (B4a 

Applying  (B2c)  yields 


1  »  B+C  +  D  +  E 


(B4b 


Applying  (B2d)  yields 


0  *  B  +  2C  +  3D  +  4E 


Applying  (B2e)  yields 


0  =  2C  +  6D  +  12E 


Finally,  (B2b)  is  applied  such  that 


U  32 (u/U  ) 

#  wr&y~  =  2C  +  6Dn  +  12En2 


which,  evaluating  at  n 


0,  yields 
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The  dimensionless  shape  parameter  is  now  defined  as 
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If  Eq  (B5)  is  substituted  into  (B4e),  the  set  of  equations  (B4a-e) 


can  be  solved  simultaneously  in  terms  of  A: 
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The  velocity  profile,  Eq  (B3),  thus  becomes 


—  =  (2n  -  2n3  +  n4)  +  ^  (n  -  3n2  + 

u  o 


3n3  -  n4)  (b 


Recall  that  the  displacement  and  momentum  thickness  are  defined  as 


Since  u  ■  Ug  for  h  >  5  ,  integration  beyond  h  =  <5  yields  nothing,  and 
the  above  equations  may  be  written 
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Uring  Eq  (B6),  the  above  integrations  can  now  be  performed,  yielding 
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The  shear  stress  at  the  wall  is  defined  as 


at  y  =  0 


Therefore,  by  using  Eq  (B6)  and  evaluating  at  n  =  0, 


T  <5 

w 


(B9) 


From  this,  the  separation  criteria,  defined  as  the  point  at  which  the 
shear  stress  at  the  wall  becomes  zero,  can  be  seen  to  be 


A  =  -12 


(BIO) 


The  last  form  of  the  momentum-integral  equation,  (Bl),  encourages 


additional  parameters  to  be  defined  as: 


Rearranging  of  terms  yields 


( B 1 6 ) 


Ue  dx  ~  l  +  fl(K  j  U~  9t  =  2  f2(K)  -  4K  -  2  K  fl(K) 

e 


One  final  parameter  is  defined  as 


F(K)  =  2  f2(K)  -  4K  -  2  K  f  (K) 


( B 1 7 ) 


By  substituting  Eq  (B17)  into  (B16)  and  then  solving  for  —  , 


£17 

g  -  |  F(K)  ♦  [4  *  fl(K)|f  5^  j  L. 
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and  solving  Eq  (B12)  for  Z  yields 
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Eqs  (B18)  and  (B19)  are  the  "working"  equations,  which  are  the  modified 
equations  that  can  be  used  to  solve  an  unsteady  boundary  layer  flow  using 


the  von  Karman-Pohlhausen  method. 


Appendix  C:  Method  of  Solution 


By  selection  of  appropriate  dimensions  and  position,  a  circle  in  one 

frame  of  reference  can  be  transformed  into  a  Joukowski  airfoil  of  desired 

thickness  and  camber  in  another  frame  of  reference  by  means  of  a  complex 

function,  commonly  known  as  a  Joukowski  transformation  (Ref  9:50-53). 

Further,  a  velocity  on  the  surface  of  a  circular  cylinder  can  be  likewise 

transformed  to  a  velocity  on  the  surface  of  the  airfoil,  such  that  for 

any  given  point  on  the  airfoil  at  any  given  angle  of  attack  a  potential 

flow  velocity  is  obtainable  (Ref  9:53-60).  By  standard  difference  methods, 

then,  the  required  partial  derivatives  of  the  potential  velocity  can  also 

be  obtained.  This  method  of  determining  geometry,  velocities,  and  partial 

derivatives,  is  described  in  depth  in  Reference  6,  and  will  not  be  repeated 

here;  however,  the  application  of  this  method  can  be  observed  in  Appendix  F. 

In  using  the  von  Karman-Pohlhausen  technique  calculation  begins  at 

9U 

the  stagnation  point.  At  this  point  it  is  known  that  =  0  and  is 

finite  and  not  equal  to  zero.  This  poses  problems  with  the  calculation  of 
clZ 

—  .  First,  the  solution  will  be  initiated  as  a  steady  flow  problem,  with 

the  pitching  motion  only  considered  after  a  finite  number  of  small  incremental 

steps,  with  negligible  effect  on  accuracy  of  the  final  solution.  Now  it 

can  be  shown  that  F(K)  must  vanish  simultaneously  with  Ue  at  the  stagnation 

dZ 

point  to  avoid  an  infinite  value  of  —  .  Solution  of  Eq  (B17)  for  F(K)  =  0 
yields  the  following  initial  values: 

A  *  7.052 

o 

K  *  0.0770 

o 
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Now  —  is  of  indeterminate  form  —  ,  which  can  be  evaluated  by  a 
dx  0  J 


limiting  process  (Ref  7:211).  The  result  is 


9u  3u 

Z  -  K  (^)  =  0.0770  (.s-S-) 

o  o  3x  o  ox  o 


and 


a2u  3u 

(— )  -  -0.0652  (r-T*)  (t-1) 

dx  o  3x  o  ox  o 


-2 


With  these  initial  values,  the  solution  can  now  be  "marched"  forward 
along  the  airfoil  as  follows: 


a.  Compute  Z,  in  the  manner  of 


Z 


1 


Z  +  (££) 
o  dx  o 


dx 


where  dx  is  an  increment  of  arc  length  along  the  airfoil. 

b.  Compute  a  new  value  of  K  from  Eq  (B 19 ) . 

c.  Compute  a  new  value  of  A  from  this  K,  either  by  solution  of 
Eq  (B12)  or  by  a  curve-fit  routine. 

d.  Compute  a  new  value  of  F(K)  and  f^(K)  by  solution  of 

Eqs  (B13),  (B14),  and  (B17).  For  the  steady  flow  solution, 
F(K)  may  be  approximated  by  the  linear  relationship 


F(K)  »  .47  -  6K  (Ref  7:213-214) 


with  reasonable  accuracy. 


>  V1  • 


e.  Compute  a  new  value  of  —  using  Eq  (B18) . 

dZ 

f.  Now,  with  new  values  for  both  Z  and  —  ,  this  procedure 
can  be  looped  back  to  step  a,  continuing  with  a  new 
incremental  step. 

The  value  of  A  (or  K)  can  be  examined  at  any  point  and  compared  with 
the  known  separation  value 


A  -  -12 


K  =>  -.1567 


For  a  given  free-stream  velocity  and  pitch  rate,  this  procedure  was 
used,  iterating  between  different  initial  angles  of  attack,  until  separation 
was  indicated  at  quarter-chord,  and  the  final  angle  of  attack  recorded. 


The  conservation  of  momentum  relationship  is  developed  by  applying 


Newton's  law  in  the  form 

F  =  ^  (mV)  =  ma 

This  analysis  is  only  valid  if  accelerations  are  measured  in  an  inertial 
reference  frame.  In  order  to  develop  a  momentum  conservation  relationship 
for  a  non-inertial  control  volume,  as  is  the  case  in  the  boundary  layer  of 
a  pitching  airfoil,  it  is  first  necessary  to  analyze  the  inertial  acceleration 
of  a  particle  in  the  control  volume  (Ref  10:109-114). 

Consider  the  problem  as  shown  in  Figure  16.  Here  the  reference  frames 
are  to  be  interpreted  as 


F  -  frame  ■  fixed  in  inertial  space 


f  -  frame  =  fixed  in  the  geometry  of  the  airfoil 


xyz  -  frame  *  fixed  in  the  control  volume 


By  the  geometry  shown  in  Figure  16,  the  coordinate  transformations  are 


=  (cos  a)f1  +  (sin  cOf^ 


F2  =  (-sin  a)f1  +  (cos  a)f2 


F3  '  {3 


=  (cos  ip)x  +  (-sin  \p)y 


f2  =  (sin  ip)x  +  (cos  ip)y 


f3  "  Z 


Combining  relationships  in  Eqs  (Dl)  and  (D2)  yields 


F^  *  cos  (a-$)  x  +  sin  (sin  ( a-ip )  y 


-sin  (a-ijO  x  +  cos  (a-^)  y 


-  " 


It  can  also  be  seen  that  the  angular  velocity  of  the  xyz  -  frame 


with  respect  to  the  F  -  frame  is 


CO 

xyz 


A 


-  a  F, 


=  -  a  z 


Now,  referring  to  Figure  17,  consider  the  motion  of  a  particle 
The  position  vector  of  the  particle  is 


P  =  R  +  r 


The  velocity  of  P  with  respect  to  the  F  -  frame  is 


dt 


Using  the  relationship  of  derivatives  of  a  vector  in  a  rotating 


(D4) 


reference  frame  (Ref  13:334-336), 
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ind  Rotating  Reference  Frame 


Differentiating  again  to  find  the  acceleration  of  P  with  respect  to 
the  F  -  frame  gives 


Fd  ,Ff\  Fd  ,*dr,  Fda)  -  -  Fdr 

dT  (  R)  +  dt  (  +  dF  *  r  *  “  *  dF 


This  reduces  to 


F —  F—  v—  —  x—  —  —  —  —  _ 

a  =  R  +  a  +  2u)  x  v  +  a)Xr  +  a)X(a)Xr) 


Newton's  law  can  now  be  applied  to  the  inertial  acceleration: 


'  f  —  —  • 

dF  =  dm  R+Xa  +  2u)X  V+u>xr  +  o)x(a}xr) 


Solving  for  acceleration  in  the  xyz  -  frame, 


dF  -  dm  R  +  20)  x  v  +  0)  x  r  +  w  x  (w  x  r)  =  dm(  a)  (D5) 


If  the  bracketed  terms  in  Eq  (D5)  are  thought  of  in  terms  of  forces, 
in  fact  body  forces  since  they  are  multiplied  by  mass,  then  this  is 
essentially  a  valid  representation  of  Newton's  law  for  a  non-inertial 
reference  frame.  Conservation  of  momentum  in  this  frame  can  now  be 
implied  in  the  usual  manner. 

The  hypothetical  body  force  in  question  then  is 


Fi 


dB  =  -  dm  R  +  2u)  x  v  +  u}xr+cox(ojxr) 


(D6) 


Referring  again  to  Figure  16,  it  can  be  seen  that,  for  a  given 
control  volume  on  the  surface  of  the  airfoil,  9  and  ip  are  constant,  the 
magnitude  of  R  is  constant,  and  a  is  considered  constant  for  the  situation 
under  investigation.  In  light  of  these  consideration,  and  given  the 
relationships  in  Eqs  (D3)  and  (D4),  each  of  the  bracketed  terms  in  Eq  (D6) 
can  be  evaluated  in  the  xyz  -  frame : 


w"  *  f  A  A  1 

R  55  R  a2  I  cos  (6+40x  -  sin  (9+^)y  J 

•  /V  A 

2u  x  *V  =  2  a  (vx  -  uy) 


(*)  x  r  =  0 


u  X  (w  X  r)  =  -  a2  (xx  +  yy) 

W * 
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Letting  <j>  =  (9  +  ip") ,  the  x-component  of  the  hypothetical  body  force 
is 

dBx  =  -  dm  |  R  a2  cos  <f>  +  2  a  v  -  a2x  J  (D7) 

which  will  be  used  in  the  analysis  of  momentum  in  the  x-direction. 

A  physical  explanation  of  the  hypothetical  body  forces  is  as  follows 

F— 

R  represents  the  component  of  acceleration  due  to  the 
translation  motion  of  the  control  volume. 

2io  x  V  represents  the  coriolis  component  of  acceleration 
due  to  the  rotational  motion  of  the  control  volume. 

(*)  x  (u)  x  r)  represe..  the  centripetal  component  of 

acceleration  due  to  the  rotational  motion 
of  the  control  volume  (Ref  13:350-353). 


Appendix  E:  Estimation  of  Mass  Introduction  Constant 


In  the  vicinity  of  the  stagnation  point,  the  frictionless  potential 
flow  velocity  distribution  is  given  by 
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as 


(Ela) 


and 
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(Ref  7:96) 


(Elb) 


where  a  is  a  constant  and  the  velocities  and  dimensions  are  as  defined  in 
Figure  18.  While  it  is  the  radial  velocity  equation  (Elb)  that  is  to  be 
used  in  the  mass  introduction  solution,  the  constant  a  can  be  approximated 
by  examination  of  the  tangential  velocity  equation  (Ela). 

Consider  the  flow  about  a  circular  cylinder  in  a  uniform  stream,  as 
depicted  in  Figure  18.  The  velocity  at  the  surface  of  the  cylinder  can  be 
shown  to  be 


u 


s 


2  U  sin  8 
00 


(Ref  8:90) 


(E2) 


For  the  region  near  the  stagnation  point,  0  is  small,  and  the 


assumption  is  made  that  sin  t  ^  8  ;  thus 
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u  -  2U  9 

s  00 


The  velocity  gradient  in  the  s-direction  is 


du  ,  du  2U 

s  _  1  s  _  00 


By  integrating  Eq  (E4)  with  respect  to  s,  the  velocity  is  obtained 
as  a  function  of  the  surface  coordinate: 


2  U  s 
00 


Comparison  of  Eqs  (Ela)  and  (E5)  yields 


for  the  given  circular  cylinder. 


The  Joukowski  airfoil  of  interest  is  obtained  by  transformation  from 
a  circular  cylinder  in  a  flow  field  as  described  in  Eq  (E2).  The  real-axis 


radius  of  this  cylinder  transforms  to  the  half-chord  of  the  Joukowski 
airfoil.  In  the  solution  method,  distances  on  the  surface  of  the  airfoil 
are  non- dimens ionali zed  with  respect  to  the  chord  of  the  airfoil.  Then 
letting  R  =  —  ,  Eq  (E6)  in  non-dimensional  form  becomes 


Appendix  F:  Computer  Program  P0HL6 


LAWRENCE.  JOHN  S. 
GAE-83D 


This  program  fs  adapted  from  program  “P0HL2* ,  Docken, 
GAE-82D,  and  will  be  used  to  continue  analysis  work  In  the 
area  of  dynamic  stall  effects.  The  representative  airfoil  1 
a  15%  Joukowskl,  pitching  at  a  constant  rate  about  Its  mid¬ 
chord  In  a  steady  freestream.  Mass  Is  being  Injected  from 
the  freestream  Into  the  boundary  layer  as  a  result  of  the 
non-Newtonian  motion  of  the  airfoil.  Symbology  within  the 
program  Is  as  follows: 

El  ■  Imaginary  unit;  !.e.,  square  root  of  -1. 

RADIUS  *  radius  of  circular  cylinder. 

AMU  *  offset  distance  of  center  of  circular  cylinder. 

ALPHA  »  angle  of  attack.  In  degrees. 

AD0T1  ■  pitch  rate,  In  degrees  per  second. 

PITCH  «  non-dimensional  pitch  rate. 

CHORD  -  chord-length  of  Joukowski  airfoil. 

ANGLE  *  radial  of  a  given  point  on  the  circular  cylinder. 
UTNF  *  freestream  velocity.  In  feet  per  second. 

CON  »  conversion  factor,  degrees  to  radians. 

X.Y.U  -  coordinates  and  potential  flow  velocity  on  airfoil 
Z,V  *  complex  coordinates  and  velocity. 

UMIN  -  ratio  of  wall  velocity  to  potential  flow  velocity. 
TIME  -  cumulative  time. 

DELT  *  Increment  of  time. 

K.N  *  Integer  counters. 

DSS  ■  Increment  of  distance  on  the  airfoil  surface. 

DUDS  *  spatial  partial  derivative  of  velocity. 

DUOT  -  partial  derivative  of  velocity  wrt  time. 

XOC  ■  chord  position. 

BETA  *  local  slope  of  airfoil  surface. 

RAS1.RAC1  ■  functions  1  of  geometry)  due  to  pitching  motion 
R.LAMDA.RK  -  Pohlhausen  shape  parameters. 

FK , ZZ , DZDS , DEL2 , F l X  ■  functions  of  shape  parameters. 

PROGRAM  P0HL6 
COMPLEX  El 

OPEN  (  lS.FILE-’FLOWIN* > 

REWIND  15 

OPEN  (  16,FILE-’FL0W0UT’ ) 

REWIND  16 
EI-(0. , 1 . ) 

RADIUS-1.131 
AMU--. 131 

READ  (15.*)  ALPHA. AD0T1 ,UINF 
ALPH1 -ALPHA 
CON-3.1415927/180. 

AD0T2-AD0T1 *CON 
THETA-180. 

T I ME*0 • 0 

CALL  DS( 180. , RAD  I  US , CON , AMU , XLE , YLE ) 

CALL  DS(0.0, RADIUS, CON, AMU , XTE , YTE ) 

XLE-ABS (XLE ) 

CHORD-XL  E  +  XTE 

P ITCH-ADOT 1 *CON*0 . 5*CHORO/UI NF 

K  - 100 
K1 -K+ 1 
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K2-(ALPHA*150)*100*K 


WRITER  1 6 , 30 )  U I  NF 
WR  I TE ( 1 6 , 40) ALP  HA 
WR!TE( 16,50 >AD0T1 
WRITEl 16.52) 

ANGLE “ALPHA ♦THETA 

CALL  U( ANGLE. RADIUS, CON, El ,UINF , AMU , ALPHA , U0  > 

CALL  DS ( ANGLE. RAD  1  US. CON, AMU, X0.V0) 

ANGLE-ANGLE-0.01 

CALL  U( ANGLE. RAD  I  US. CON, El . U I NF , AMU . ALPHA . U 1 ) 

CALL  DS! ANGLE, RAD  I  US, CON, AMU. XI ,Y1 ) 

ANGLE-ANGLE-0.01 

CALL  U( ANGLE. RAD  I  US. CON, El , U I NF , AMU , AL P HA , U2 ) 

CALL  OS ! ANGLE, RADI  US, CON, AMU. X2.Y2) 

DS2- ( SORT  < (X2-X1 ) **2  +  ! Y2-Y1 >--2 > > /CHORD 
DS 1 - ( SORT ( (X1-X0>*"2*1Y1-Y0>**2 ) ) /CHORD 

Stagnation  point  velocity  gradient  computed  using  a 
forward  difference  method:  all  other  velocity  gradients 
computed  using  central  difference  method. 

DUDS-(U2-U0)/(DS1*DS2) 

Second  derivative  of  velocity  computed  using  a 
Taylor’s  series  expansion. 

D2UDS2-<U2-2.-Ul"-U0)/<  <DSl-*DS2>/2.  )-»2 

Enter  initial  boundary  layer  parameters. 

RLAMDA-7 .052 

RK-0.0770 

FK-0.0 

DZDS--0.0652*D2UDS2/<DUDS**2> 

ZZ-RK/DUDS 

N-50 

ANGLE-ALPHA+THETA-0.01 
XOC- ( X0+XLE  > /CHORD 

WRITEl 16,1) XOC , U0 , OUDS , D2UDS2 ,FK.RK,ZZ,DZDS 

ADOT -0.0 
DO  10  J-l  ,K 

Function  of  this  loop  Is  to  compute  boundary  layer 
parameters  at  stagnation  point,  allowing  the 
boundary  layer  to  steady-out  before  subjecting  It 
to  a  pitching  motion. 

N-N*  1 


Compute  pertinent  boundary  layer  parameters. 

ZZ-DZDS-DS1+ZZ 
RK-ZZ*DUDS 
FK- . 47-6 . *RK 
DZDS-FK/U 1 

DELT-CHORD*DS 1 / U 1 
TIME-TIME+DELT 

CALL  U (ANGLE, RAD  I  US, CON, El .UINF, AMU. ALPHA. U2> 
DUDT-IU2-U1  l/DELT 
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ANGLE-ANGLE-0.01 
ANGLE  1 -ANGLE-0.01 

CALL  U1ANGLE1 , RAD  I  US . CON . E I . U I NF . AMU , AL PHA . U2 > 

CALL  DS ( ANGLE  1 . RAO  I  US . CON , AMU . X2 . Y2 ) 

ANGLE0- ANGLE +0.01 

CALL  U ( ANGLE0 . RAD  I  US, CON, El ,  U  I  NF .AMU .ALPHA, U0> 

CALL  OS ! ANGLE0 .RAO  I  US , CON .AMU , X0 , Y0 ) 

CALL  Ut ANGLE, RADIUS. CON. El ,UINF .AMU .ALPHA . U 1 ) 

CALL  OS (ANGLE, RAD  I  US, CON. AMU, XI ,Y1  ) 

DS 1 - ( SORT ( (X1-X0)**2+(Y1-Y0)**2) ) /CHORD 
DS2-(SQRT( (X2-X1 >**2+< Y2-Y1 >**2> ) /CHORD 

Compute  the  arc  length  and  the  velocity  gradient. 

DSS-DS 1+DS2 
DUDS-(U2-U0>/DSS 
XOC-(Xl+XLE)/CHORD 
IF1N.LT. 50)  GO  TO  10 

N-0 

WRITE1 16, 1 )XOC,Ul ,DUDS.DUDT.FK.RK,ZZ,DZDS 
CONTINUE 

ADOT-ADOT1 

N-0 

RAS 1-0.0 
RAC  1-0.0 
DO  20  O-Kl ,K2 

Function  of  this  loop  Is  to  compute  the  behavior 
of  the  boundary  layer  as  It  Is  subjected  to  a 
pitching  motion. 

N-N+l 

Compute  the  pertinent  boundary  layer  parameters. 

ZZ-DZDS*0S 1 +ZZ 

RK-ZZ* (  DUDS* ( 1 .+RAS1/U1 >*DUDT/U1 ) 

CALL  POHL  <  RK , RLAMOA ) 

DEL2-37./31 5. -RLAMDA/945 . -1RLAMDA--Z 1/9072. 

FK*2.*DEL2M2.-. 3683*RLAMDA+ .0104-RLAMDA--2* 
<RLAMDA*-3>/4536. ) 

F 1 K- ( . 3-RLAMDA/ 1 20. 1/DEL2 

DZDS-(FK+(4.+FlK)*ZZ*DUDT/Ul+(4.+2.*FlK-2./DEL2)*ZZ 
•RAS1 *DUDS/U1-8.*ZZ*RAC1/DEL2 l/Ul 

Compute  the  time  Increment  for  a  particle  to 
travel  from  point  (1)  to  point  11+1). 

DELT-CH0RD*DS1/U1 

TIME-T1ME+DELT 

DALPHA-DELT*ADOT 

ANGLE-ANGLE+DALPHA 

ALPH1-ALPH1+DALPHA 

CALL  U( ANGLE, RAD  I  US, CON, El , U I NF , AMU , ALPH 1 ,U2C> 

CALL  N I U (ANGLE . RAD  I  US , CON , AMU , AD0T2 , X 1 , Y 1 . U2C , U2 , RAS2 .RAC 2 ) 

Compute  the  unsteady  velocity  gradient. 

0UDT-(U2-U1 >/DELT 
ANGLE-ANGLE-0.01 
ANGLE  1 -ANGLE-0.01 

CALL  U< ANGLE  1 .RADI US, CON. El .UINF .AMU, ALPH 1 ,U2C> 

CALL  0S(ANGLE1 , RAD  I  US . CON , AMU , X2 . Y2 ) 


CALL  N I U ! ANGLE  1 .RADIUS, CON .AMU ,ADOT2,X2,Y2,U2C,U2.RAS2. RAC 2 ) 
ANGLE0-ANGLE+0.01 

CALL  U ( ANGLE#, RAD  I  US . CON , E I .  UINF , AMU. AL PHI , U0C > 

CAL L  OS  <  ANGLE# . RAO  I  US , CON . AMU , X0 , Y0  > 

CALL  N I U ( ANGLE# . RAD  I  US , CON .AMU . AOOT2 , X# , Y0 , U0C , U0 , RAS# , RAC#  > 
CALL  U< ANGLE. RAD  I  US, CON. El .UINF . AMU , ALPH 1 . U 1 C > 

CALL  OS (ANGLE. RAO  I  US. CON. AMU, XI ,Y1 ) 

CALL  NIU(ANGLE.RADIUS.CON,AMU.ADOT2,Xl ,Y1 .U1C.U1 .RAS1 ,RAC1 ) 

Compute  arc  length  and  velocity  gradient. 

DS 1  “  ( SORT ( !X1-X0)*“2+!Y1-Y0)**2) ) /CHORD 

DS2-(SQRT( (X2-X1 > *»2* < Y2-Y1 > *»2 ) ) /CHORD 

DSS*DS 1 +DS2 

DUDS-( U2-U#>/DSS 

UDU0S-U1*0U0S 

XOC*!Xl+XLE) /CHORD 

Stop  the  computation  at  the  quarter -chord . 

IF(XOC.GE. 0.250)  GO  TO  25 
IF( N. LT. 250)  GO  TO  20 

N*0 

WRITE1 16, 1 )XOC.Ul .DUDS.DUDT.FK.RK.ZZ.DZDS 
20  CONTINUE 

25  WRITE! 16, 1 1XOC.U1 .DUDS.DUDT.FK.RK.ZZ.DZDS 

WRITE1 16.45  > ALP  HI 
WRITE ( 16.55) PITCH 
WRITEl 16,60)RK 
WRITE1 16,80) 

WRITEl 16,81 )TIME 
UMIN-RAS1/U1 
WRITEl 16 ,85  >UM IN 

1  FORMAT!4X,F6.3.2(4X,F10.3> ,4X,F7.3.4X,F7.4.4X.F8.4.2(  4X.E9.3) > 

30  FORMAT1 1H1 , ’BOUNDARY-LAYER  PARAMETERS  FOR  ", F6 . 2 ." FT/SEC "/ > 

40  FORMAT ( "  INITIAL  ANGLE  OF  ATTACK!  ’.F6.3.*  DEGREES"/) 

45  FORMAT ( / *  FINAL  ANGLE  OF  ATTACK:  ".F6.3."  DEGREES"/) 

50  FORMAT ( "  PITCH  RATE:  ’.F7.3,"  DEGREES/SEC"/) 

52  FORMAT ( 6X , "XOC " • 10X, "U" , 1 IX  ,  " DUDS " , 9X , "DUDT" , 8X . " FK" . 

♦  9X,"RK",11X,"Z’,11X,"DZDS’/) 

55  FORMAT! "  PITCH  PARAMETER:  ",F7.5/> 

60  FORMAT ( "  K  AT  THE  QUARTER-CHORD:  “,F8.4/> 

80  FORMAT!"  TIME  TO  REACH  THE  QUARTER-") 

81  FORMAT!"  CHORD  FROM  THE  STAGNATION  POINT:  ",F7.5,"  SEC"/) 

85  FORMAT!"  UWALL/UE  AT  THE  QUARTER-CHORD:  ",E9.3/> 

STOP 

100  END 


SUBROUTINE  U ! ANGLE ,RAD I  US , CON , E I . U I NF , AMU .ALPHA . UU > 
COMPLEX  CMPLX.Z.EI .DZETA.W 

Function  of  this  subroutine  Is  to  to  compute  the  local 
value  of  velocity  on  a  Joukowsk <  airfoil  using  complex 
potential  flow  theory. 

X«RADIUS*COS!ANGLE*CON) 

Y-RADIUS*SIN!ANGLE*CON) 

Z-CMPLX1X.Y) 

W»U INF*! ( 1 . ,0. ) - ! RAD IUS**2)/Z**2» 

♦  ! 2 . *E I *RAD I US*S I N  <  ALPHA*CON  )  >/Z> 

X-X+AMU 
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* 

* 


Z  changed  to  represent  values  of  coordinates  used  In 
the  transf ormat 1  on  equation. 


Z"CMPLX ( X  V) 

DZETA"(Z**2-<  RAO  1 US*AMU )*"2)/Z**2 
UU-CABS ( W ) /CABS ( DZETA) 

RETURN 

END 


SUBROUTINE  DS ! ANGLE . RAD  I  US . CON . AMU . X . Y > 

COMPLEX  CMPLX.Z 

Function  of  this  subroutine  is  to  compute  the  arc 
length  between  two  points  on  the  Uoukowskl  airfoil. 

X-RAD I US»COS  <  ANGL  E  *CON ) *AMU 
Y«RADIUS*SIN<ANGLE*CON> 

Z-CMPLX(X , Y> 

Z-Z-M (RADI US* AMU ) **2  )  /Z 
X-REAL(Z) 

Y-AIMAG(Z) 

RETURN 

END 


SUBROUTINE  POHL(RK.RLAMDA) 

Function  of  this  subroutine  Is  to  compute  the  value 
of  the  separation  parameter,  lamda,  given  a  value 
of  K,  as  computed  In  the  main  program. 


RK1  — . 

160 

RK2— . 

112 

RK3"0. 

00 

RK4-0. 

06 

RK5"0 . 

076 

RK6-0. 

086 

RK7-0. 

0949 

1F(RK. 

LE.RK1  ) 

GO 

TO 

10 

IF (RK. 

LE.RK2) 

GO 

TO 

20 

I F  ( RK  . 

LE.RK3) 

GO 

TO 

30 

IF  IRK. 

LE.RK4) 

GO 

TO 

40 

I F  (  RK . 

LE.RK5) 

GO 

TO 

50 

IF1RK. 

LE.RK6) 

GO 

TO 

60 

I F  ( RK  . 

GT.RK7) 

GO 

TO 

70 

RLAMDA 

t"  .01  49*’ 

'2- 

(RK-0.1 

R LAMDA* 1 2 . - 100. "SORT ( RLAMDA ) 
RETURN 

10  RLAM0A-(2./.Br:2)*RK>  W.0 

RETURN 

20  RLAMDA" (4./.044)"RK+2. 18 
RETURN 

30  RLAMDA" ( 10./.14>*RK 

RETURN 

40  RLAMDA- 8 3 . 33*RK 

RETURN 
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J-J- 


RLAMOA-- 1. 9+115. *RK 
RETURN 


50 


60  RLAMDA--6. 54+176. *RK 
RETURN 

70  RLAMDA-12. 

RETURN 

END 


SUBROUTINE  NIU{ANGLEC,RADIUS.CON,AMU.ADOT2,XC,YC,UC,UNI  . 

+  RASF.RACFl 

Function  of  this  subroutine  Is  to  compute  the  value  of  U  In 
the  non-lnertlal  reference  frame,  and  the  values  of  the 
required  functions  of  geometry. 


* 


R-SORT<XC*«2+YC*«2> 
BETl-ATAN(YC/(-XC) ) 

ANGLEM-ANGLEC+ . 1 
ANGLEP-ANGLEC-  .  1 

CALL  DS ( ANGLEM . RAD  I  US , CON .AMU , XM . YM ) 
CALL  DS (ANGLER .RADIUS .CON , AMU .XP ,YP 1 
DY-ABS( YP-YMl 
DX-ABS(XP-XM) 

Q«1000.*DX 

IF ( DY.GT.Q1THEN 

BETA-1.5707963 

ELSE 

BETA-BET1 +ATAN ( DY/DX ) 

END  IF 

RASF -R*ADOT2*S I N ( BETA  > 

RACF -R"ADOT2*COS ( BETA ) 

UNI-UC-RASF 

RETURN 

END 
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extended  to  the  second  case.  Using  the  Moore-Rott-Sears  model  for  flow 
separation  criteria,  analytical  results  were  computed  and  compared  with 
experimental  data.  Reduction  in  adversity  of  the  pressur'  gradient 
accounts  for  only  a  fraction  of  the  total  dynamic  effect,  rnd  it  is 
proposed  that  mass  introduction  into  the  boundary  layer  from  the 
free-stream  may  be  a  strongly  contributing  factor.  This  phenomena  is 
demonstrated  to  have  a  large  effect,  and  an  argument  is  presented  for 
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